Catch that asteroid!¶

First, we need to increase the timeout time to allow the download of data occur properly

[1]:

from astropy.utils.data import conf
conf.dataurl

[1]:

'http://data.astropy.org/'

[2]:

conf.remote_timeout

[2]:

10.0

[3]:

conf.remote_timeout = 10000

Then, we do the rest of the imports and create our initial orbits.

[4]:

from astropy import units as u
from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris
solar_system_ephemeris.set("jpl")

from poliastro.bodies import *
from poliastro.twobody import Orbit
from poliastro.plotting import StaticOrbitPlotter
from poliastro.plotting.misc import plot_solar_system

EPOCH = Time("2017-09-01 12:05:50", scale="tdb")

WARNING: AstropyDeprecationWarning: astropy.extern.six will be removed in 4.0, use the six module directly if it is still needed [astropy.extern.six]

[5]:

earth = Orbit.from_body_ephem(Earth, EPOCH)
earth

[5]:

1 x 1 AU x 23.4 deg (ICRS) orbit around Sun (☉) at epoch 2017-09-01 12:05:50.000 (TDB)

[6]:

earth.plot(label=Earth);

/home/lobo/Github/poliastro/src/poliastro/twobody/orbit.py:1163: UserWarning:

Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned

../_images/examples_Catch_that_asteroid!_8_1.png

[7]:

florence = Orbit.from_sbdb("Florence")
florence

[7]:

1 x 3 AU x 22.1 deg (HeliocentricEclipticIAU76) orbit around Sun (☉) at epoch 2458600.5008007586 (TDB)

Two problems: the epoch is not the one we desire, and the inclination is with respect to the ecliptic!

[8]:

florence.rv()

[8]:

(<Quantity [-2.76132873e+08, -1.71570015e+08, -1.09377634e+08] km>,
 <Quantity [13.17478676, -9.82584124, -1.48126637] km / s>)

[9]:

florence.epoch

[9]:

<Time object: scale='tdb' format='jd' value=2458600.5008007586>

[10]:

florence.epoch.iso

[10]:

'2019-04-27 00:01:09.186'

[11]:

florence.inc

[11]:

$22.142394 \; \mathrm{{}^{\circ}}$

We first propagate:

[12]:

florence = florence.propagate(EPOCH)
florence.epoch.tdb.iso

[12]:

'2017-09-01 12:05:50.000'

And now we have to convert to the same frame that the planetary ephemerides are using to make consistent comparisons, which is ICRS:

[13]:

florence_icrs = florence.to_icrs()
florence_icrs.rv()

[13]:

(<Quantity [ 1.46404253e+08, -5.35752831e+07, -2.05656912e+07] km>,
 <Quantity [ 7.34329037, 23.47561546, 24.12063695] km / s>)

Let us compute the distance between Florence and the Earth:

[14]:

from poliastro.util import norm

[15]:

norm(florence_icrs.r - earth.r) - Earth.R

[15]:

$6967159.9 \; \mathrm{km}$

This value is consistent with what ESA says! $7\,060\,160$ km

[16]:

abs(((norm(florence_icrs.r - earth.r) - Earth.R) - 7060160 * u.km) / (7060160 * u.km))

[16]:

$0.013172521 \; \mathrm{}$

[17]:

from IPython.display import HTML

HTML(
"""<blockquote class="twitter-tweet" data-lang="en"><p lang="es" dir="ltr">La <a href="https://twitter.com/esa_es">@esa_es</a> ha preparado un resumen del asteroide <a href="https://twitter.com/hashtag/Florence?src=hash">#Florence</a> 😍 <a href="https://t.co/Sk1lb7Kz0j">pic.twitter.com/Sk1lb7Kz0j</a></p>&mdash; AeroPython (@AeroPython) <a href="https://twitter.com/AeroPython/status/903197147914543105">August 31, 2017</a></blockquote>
<script src="//platform.twitter.com/widgets.js" charset="utf-8"></script>"""
)

[17]:

La @esa_es ha preparado un resumen del asteroide #Florence 😍 pic.twitter.com/Sk1lb7Kz0j
— AeroPython (@AeroPython) August 31, 2017

And now we can plot!

[18]:

frame = plot_solar_system(outer=False, epoch=EPOCH)
frame.plot(florence_icrs, label="Florence");

/home/lobo/Github/poliastro/src/poliastro/twobody/orbit.py:1163: UserWarning:

Frame <class 'astropy.coordinates.builtin_frames.icrs.ICRS'> does not support 'obstime', time values were not returned

../_images/examples_Catch_that_asteroid!_26_1.png

The difference between doing it well and doing it wrong is clearly visible:

[19]:

frame = StaticOrbitPlotter()

frame.plot(earth, label="Earth")

frame.plot(florence, label="Florence (Ecliptic)")
frame.plot(florence_icrs, label="Florence (ICRS)");

../_images/examples_Catch_that_asteroid!_28_0.png

We can express Florence’s orbit as viewed from Earth. In order to do that, we must set the Earth as the new attractor by making use of the change_attractor() method. However Florence is out of Earth’s SOI, meaning that changing the attractor from Sun to Earth has no physical sense. We will make use of force=True argument so this method runs even if we know that we are out of new attractor’s SOI.

[20]:

florence_hyper = florence.change_attractor(Earth, force=True)

/home/lobo/Github/poliastro/src/poliastro/twobody/orbit.py:489: PatchedConicsWarning:

Leaving the SOI of the current attractor

Previous warning was raised since Florence’s orbit as seen from Earth is hyperbolic. Therefore if user wants to propagate this orbit along time, there will be some point at which the asteroid is out of Earth’s influence (if not already).

We now retrieve the ephemerides of the Moon, which are given directly in GCRS:

[21]:

moon = Orbit.from_body_ephem(Moon, EPOCH)
moon

[21]:

367937 x 405209 km x 19.4 deg (GCRS) orbit around Earth (♁) at epoch 2017-09-01 12:05 (TDB)

[22]:

moon.plot(label=Moon);

../_images/examples_Catch_that_asteroid!_34_0.png

And now for the final plot:

[23]:

import matplotlib.pyplot as plt

frame = StaticOrbitPlotter()

# This first plot sets the frame
frame.plot(florence_hyper, label="Florence")

# And then we add the Moon
frame.plot(moon, label=Moon)

plt.xlim(-1000000, 8000000)
plt.ylim(-5000000, 5000000)

plt.gcf().autofmt_xdate()

/home/lobo/Github/poliastro/src/poliastro/twobody/orbit.py:1104: OrbitSamplingWarning:

anomaly outside range, clipping

../_images/examples_Catch_that_asteroid!_36_1.png

Per Python ad astra!